COMMENT ⊗   VALID 00002 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00002 00002	\magnification\magstephalf
C00011 ENDMK
C⊗;
\magnification\magstephalf
\input macro.tex
\def\today{\ifcase\month\or
  January\or February\or March\or April\or May\or June\or
  July\or August\or September\or October\or November\or December\fi
  \space\number\day, \number\year}
\baselineskip 14pt
\rm
\line{\sevenrm a16.tex[106,phy] \today\hfill}

\bigskip
\ctrline{[Subject: Precision.]}

\bigskip
\ctrline{\bf Testing a Numerical Algorithm}

\bigskip
Programmers rely on libraries of subalgorithms to perform standard computations.
Some are built into the translators for programming languags, like the
standard functions of Pascal. Others are available on text files, written
in common programming languages, and may be copied into programs. A~typical
computer installation might have standard algorithms to solve simultaneous
linear equations, integrate numerically, convert between polar and
rectangular coordinates, perform arithmetic on complex numbers, and many 
others. Where precision is important, it is unwise to use such algorithms
without investigating their precision.

If the desired results of a subalgorithm are characterized by a property
that can be tested, a~program can usually evaluate the error in the results.
Look at the standard function SQRT$(x)$ for computing approximate square
roots. The true square root, $y=\sqrt{x}\,$, is characterized by
$y↑2=x$, $y≥0$, which is readily tested. We can not expect this test
to be satisfied exactly by SQRT$(x)$, though, because of rounding
error. We expect SQRT$(x)=\sqrt{x}\,(1+e)$, with $|e|$ small, and a program
can use the characterization of the true square root to estimate~$e$.

If we could do exact arithmetic, we could compute
$${\bigl(\hbox{SQRT}(x)\bigr)↑2/x-1\over 2}={(1+e)↑2-1\over 2}
={2e+e↑2\over 2}=e+{e↑2\over 2}\,,$$
which is almost exactly $e$.

\def\boxbinary#1%
{\hbox to 1em{\hbox to 1em{\hfil#1\hfil}\hskip-1.05em minus1em
\vbox to 9.8222pt{\hrule height .5pt
\hbox to 1.2em{\vrule height 1.2em width .5pt\hfil\vrule height 
1.2em width .5pt}
\hrule height .5pt\vss}\hskip-.1emminus1em}}

\def\boxbinarya#1%
{\hbox to .5em{\hbox to .5em{\hfil#1\hfil}\hskip-.525em minus.5em
\vbox to 4.9111pt{\hrule height .25pt
\hbox to .6em{\vrule height .6em width .25pt\hfil\vrule height 
.6em width .25pt}
\hrule height .25pt\vss}\hskip-.05emminus.5em}}

\def\boxp{\mathbin{\boxbinary{$+$}}}
\def\boxm{\mathbin{\boxbinary{$-$}}}
\def\boxd{\mathbin{\boxbinary{$/$}}}
\def\boxa{\mathbin{\boxbinary{$\ast$}}}
\def\boxtwo{\boxbinarya{$\scriptstyle 2$}}

The next best thing is to do a similar computation with REAL arithmetic,
and take into account the errors introduced by the testing itself.
If we write
$$\eqalign{y&:=\hbox{SQRT}(x)\cr
D&:=\bigl(\hbox{SQR}\bigl(\hbox{SQRT}(x)\bigr)-x\bigr)/(2\,\ast\,x)\cr}$$
we actually get
$$\eqalign{y&:=\sqrt{x}\,(1+e)\cr
D&:=\qquad\qquad\bigl(\bigl(\sqrt{x}\,(1+e)\bigr)↑{\boxtwo}\,\boxm x\bigr)
\;\boxd\;(2\,\boxa\,x)\cr
\noalign{\medskip}
&={\bigl(\bigl(\sqrt{x}(1+e)\bigr)↑2(1+e↓1)-x\bigr)(1+e↓2)\over 2x(1+e↓3)}(1+e↓4)
\cr}$$
where, on a binary machine, $2x$ will be computed exactly, so $e↓3$ is~0, and
$$\eqalign{D&={\bigl(x(1+e)↑2(1+e↓1)-x\bigr)\over 2x}(1+e↓2)(1+e↓4)\cr
&=e+{e↓1\over 2}+\hbox{negligible terms like }{e↓1e↓2\over 2}\,.\cr}$$
Since $|e↓1|≤\hbox{EPSILON}/2$ (assuming a rounding machine), we find $e$
must lie between $D-\hbox{EPSILON}/4$ and $D+\hbox{EPSILON}/4$.

Subtracting EPSILON$/4$ from the absolute value of $D$, we can print
$\hbox{ABS}(D)-\hbox{EPSILON}/4$ in the knowledge that it is a good estimate
of the relative error in SQRT, that it never overestimates the error, but
that it may underestimate the error by $\hbox{EPSILON}/4$. If we try this
on a substantial sample of values of~$x$, and the printed results are
never more than $\hbox{EPSILON}/2$ (or at most EPSILON), we have some reason
to think SQRT is about as precise as we can hope for.

\smallskip
\noindent{\bf Exercise:} On a decimal machine, when is $2\,\ast\,x$ computed
exactly? 

\smallskip
\noindent{\bf Answer:} When the mantissa is less than 5.0, or when its last
digit is 0 or~5.

\smallskip
\noindent{\bf Exercise:} How does the analysis of SQRT change on a truncating
machine? 

\smallskip
\noindent{\bf Answer:} Since $-\hbox{EPSILON}<e↓1≤0$, we find $e$ must lie
between $D$ and $D+\hbox{EPSILON}/2$.

\bigskip
\ctrline{\bf Analysis of error in Vancouver stack index.}

\smallskip
At each stock transaction, a correction to the index, $C↓i$, is computed,
truncated to three decimal places, and added to the index. Even if the
computation of $C↓i$ and the addition were both exact, the truncation would
subtract a number between 0 and~0.001, on the average 0.0005. After $N$
corrections, the computed result is below the true result by the sum of $N$
such numbers, or about $0.0005N$. To reach an error around 500~points,
$N$~must be about $500/0.0005=10↑6$. To reach this level in 1.5~years
$\approx 500$~days would require $10↑6/500=2000$ transactions per day, on the
average. If this number is too high, there may be other errors in the program,
and we may be hearing of more trouble in Vancouver. Keep reading the {\sl Wall
Street Journal\/} for the latest news on numerical precision.

\vfill\eject\end